!> 求解域
module sph_region

    use sph_kinds, only: rk
    use sph_particle, only: particle_t
    use sph_rigid_body, only: rigid_body
    use sph_env, only: dim
    use sph_terminal

    type region_t
        integer :: Nreal, Nvirt, Nrigid_particle, Ntotal, Nrigid  !! 粒子数
        real(rk) :: hsml  !! 平滑长度
        type(particle_t) :: particles  !! 全体粒子
        type(particle_t) :: virtual_particles  !! 虚拟粒子，用于恢复粒子位置
        type(rigid_body), allocatable :: rigid(:)  !! 刚体
    contains
        procedure :: get_number_of_rigid_body
        procedure :: load_pif_nml
        procedure, private :: check_dim
    end type region_t
    character(*), parameter :: pif_nml = 'pif.nml'
    character(*), parameter :: pif_h5part = 'pif.h5part'
contains

    !> 获取刚体数目
    pure function get_number_of_rigid_body(this) result(nrigid)
        class(region_t), intent(in) :: this
        integer :: nrigid
        nrigid = this%Nrigid
    end function get_number_of_rigid_body

    !> 保存虚拟域粒子中可能会被改变的量
    !> @note 目前仅保存位置、速度、密度
    pure subroutine save_virtual_particles(real_part, virtual_part, nvirt, nreal)
        type(particle_t), intent(in) :: real_part  !! 真实粒子
        type(particle_t), intent(inout) :: virtual_part  !! 虚拟粒子
        integer, intent(in) :: nvirt, nreal  !! 粒子数

        virtual_part%loc(:, 1:nvirt) = real_part%loc(:, nreal + 1:nreal + nvirt)
        virtual_part%vel(:, 1:nvirt) = real_part%vel(:, nreal + 1:nreal + nvirt)
        virtual_part%rho(1:nvirt) = real_part%rho(nreal + 1:nreal + nvirt)
        virtual_part%p(1:nvirt) = real_part%p(nreal + 1:nreal + nvirt)

    end subroutine save_virtual_particles

    !> 更新虚拟域粒子中可能会被改变的量
    !> @note 目前仅更新位置、速度、密度。
    pure subroutine update_virtual_particles(real_part, virtual_part, nvirt, nreal)
        type(particle_t), intent(inout) :: real_part  !! 真实粒子
        type(particle_t), intent(in) :: virtual_part  !! 虚拟粒子
        integer, intent(in) :: nvirt, nreal  !! 粒子数

        real_part%loc(:, nreal + 1:nreal + nvirt) = virtual_part%loc(:, 1:nvirt)
        real_part%vel(:, nreal + 1:nreal + nvirt) = virtual_part%vel(:, 1:nvirt)
        real_part%rho(nreal + 1:nreal + nvirt) = virtual_part%rho(1:nvirt)
        real_part%p(nreal + 1:nreal + nvirt) = virtual_part%p(1:nvirt)

    end subroutine update_virtual_particles

    !> 从pif文件中读取参数
    subroutine load_pif_nml(this, nml_file, nums, hsml)
        class(region_t), intent(inout) :: this
        character(*), intent(in) :: nml_file
        integer :: nums(5)  !! 整数变量
        real(rk) :: hsml  !! 光滑长度
        integer :: iunit, real_kind, dim, nrigid, index, ik
        allocatable :: index(:, :)
        character(6) :: creator
        character(23) :: create_date
        namelist /pif_attr/ creator, create_date, real_kind, nums, dim, hsml
        namelist /rigid/ nrigid, index
        open (newunit=iunit, file=nml_file, action='read', status='old')
        read (iunit, nml=pif_attr)
        this%Nreal = nums(1)
        this%Nvirt = nums(2)
        this%Nrigid_particle = nums(3)
        this%Ntotal = nums(4)
        this%Nrigid = nums(5)
        this%hsml = hsml
        call this%check_dim(dim)
        if (this%Nrigid > 0) then
            allocate(this%rigid(this%Nrigid), index(2, this%Nrigid))
            read (iunit, nml=rigid)
            do ik = 1, this%Nrigid
                this%rigid(ik)%istart = index(1, ik)
                this%rigid(ik)%iend = index(2, ik)
            end do
        end if
        close (iunit)
    end subroutine load_pif_nml
    !> 检查维度
    subroutine check_dim(this, dimension)
        class(region_t), intent(in) :: this
        integer, intent(in) :: dimension
        if (dim /= dimension) then
            call term%error('Dimension mismatch!')
        end if
    end subroutine check_dim
end module sph_region
